Nuclear quantum effects in ab initio dynamics: 
theory and experiments for lithium imide 
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Owing to their small mass, hydrogen atoms exhibit strong quantum behavior even at room tem- 
perature. Including these effects in first principles calculations is challenging, because of the huge 
computational effort required by conventional techniques. Here we present the first ab-initio ap- 
plication of a recently- developed stochastic scheme, which allows to approximate nuclear quantum 
effects inexpensively. The proton momentum distribution of lithium imide, a material of interest 
for hydrogen storage, was experimentally measured by inelastic neutron scattering experiments and 
compared with the outcome of quantum thermostatted ab initio dynamics. We obtain favorable 
agreement between theory and experiments for this purely quantum mechanical property, thereby 
demonstrating that it is possible to improve the modelling of complex hydrogen-containing materials 
without additional computational effort. 



Nuclear quantum effects play an important role in de- 
termining the properties of compounds containing light 
elements, hydrogen in particular. In order to assist the 
interpretation of experiments, accurate theoretical mod- 
eling is highly desirable. Unfortunately, conventional 
techniques [TJ [2] can be orders of magnitude more expen- 
sive than methods which treat the nuclei as classical par- 
ticles. As a consequence, in ab initio simulations nuclear 
quantum effects have seldom been included |3j 0]. 

A stochastic molecular dynamics framework based on 
generalized Langevin equations has been recently de- 
vised. Among the many possible applications [5 -7 , it al- 
lows one to model to a good approximation nuclear quan- 
tum effects at negligible additional effort with respect 
to purely classical dynamics. Preliminary tests based 
on empirical force fields demonstrated satisfactory agree- 
ment with path integral and experimental results [6[ [7]. 

An additional advantage of this approach is that not 
only atomic configurations, but also the momentum re- 
produce the quantum distribution. On the contrary, 
computing the momentum distribution within a path in- 
tegral formalism [2] [S] requires special techniques and, de- 
spite recent developments [9], further increases the com- 
putational effort. In the classical limit, the distribution 
of the momentum p of a particle is Gaussian, n(p) oc 
exp(— p 2 /2m/c#T), and depends only on the temperature 
T and the particle's mass m. Conversely, in a quantum 
mechanical description n(p) reflects the local potential 
experienced by the particle. Deviation of the proton mo- 
mentum distribution (PMD) from the classical one is a 
very sensitive probe of the quantum-mechanical behav- 
ior of hydrogen atoms. Experimental measurements of 
the PMD have been made feasible since the advent of 
spallation neutron sources. Indeed, the intense fluxes of 
neutrons in the 1-100 eV energy range provided by these 



facilities allows one to study the short time (10 -16 s) 
dynamical properties of the proton in different hydro- 
gen containing systems, as well as quantum fluids and 
solids [10]. 

In this Letter we apply for the first time the "quantum 
thermostat" together with ab initio molecular dynamics, 
in order to model lithium imide. Besides its technologi- 
cal significance as a material for hydrogen storage [TTJ [12], 
Li2NH is well suited as a benchmark. In fact, the pres- 
ence of libration modes of the NH bonds is likely to in- 
troduce significant anharmonicities, which rule out the 
possibility of an accurate treatment by harmonic lattice 
dynamics. 

Theoretical results are compared with the 
experimentally-determined PMD in lithium imide, 
which was obtained by means of Deep Inelastic Neutron 
Scattering (DINS) measurements on the VESU- 
VIO spectrometer at the ISIS spallation neutron source 
(Rutherford Appleton Laboratory, United Kingdom) [13J, 
in the Resonance Detector (RD) configuration and using 
the Foil Cycling Technique [14] [15] that provides a nar- 
row resolution suitable for line shape analysis on PMD. 
The high energy and wave vector transfers achievable 
with DINS allow one to describe the scattering event 
within the framework of the impulse approximation 
(I A) with a very good degree of accuracy [TO] ITS] , so as 
to extract the proton momentum distribution directly 
from the experimental data. Actually the value of the 
wave vector transfer at the maximum of the proton 
recoil peak ranges from 35 A" 1 to 150 A" 1 for the 
complete set of detectors used in the present experiment. 
Following Refs.fTT] ITS] , one can easily evaluate the 
coefficient which multiplies the first term beyond the 
impulse approximation, namely the third derivative 
of the IA response function itself. Using appropriate 
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FIG. 1. A representation of the initial configuration of the 
atoms in the supercell used for our simulations [25]. The 
stable, tetrahedral clusters of interstitial Li atoms are high- 
lighted. 




physical quantities for lithium imide, one finds (in 
the low temperature limit and for the aforementioned 
values of the wave vector transfer) that this "final state 
effect" coefficient lies in between 3.67 A -3 and 0.856 
A -3 . This ensures that the impulse approximation is 
already reached during the reported neutron scattering 
measurement to any practical purpose. 

The polycrystalline sample of I^NH was synthesized 
by thermal decomposition of commercial lithium amide 
( Sigma- Aldrich Inc., reagent grade) at T = 623 K and 
p = 10 -3 Pa for 4 h, according to the reaction LiNH2 — >• 
\ Li2NH + \ NH3 [IHEO], in a furnace equipped with 
turbomolecular vacuum pump. X-ray powder diffraction 
measurements (CuKq, radiation) showed the sample to be 
well crystallized and to contain a small quantity of I^O, 
which was already present in the pattern of the starting 
lithium amide. This minor contamination was reported 
also in the previous studies [T9] [20] ; on the other hand, 
no traces of LiOH or of other impurities were observed. 

Being a powdered sample, DINS measures the spher- 
ically averaged PMD. The detailed experimental proce- 
dures to extract the PMD from DINS data can be found 
in Refs.[21-24 where data reduction and line shape anal- 
ysis procedures are described in details. 

Simulations were performed using a supercell contain- 
ing 192-atoms (Figure [T]), which were arranged according 
to the partially-disordered structure which was recently 
proposed as a model of the low-temperature, Fd3m phase 
of Li 2 NH[25 . In this structure, Li atoms occupy the 
tetrahedral sites of the fee lattice of nitrogen atoms. Two 
types of Li vacancies with different local symmetry are 
present. One kind of vacancy is tetrahedrally coordinated 
by N-H groups, while the second one is tetrahedrally 
coordinated by Li interstitials. These tetrahedral clus- 
ters of Li interstitials are found to stabilize the structure 
considerably, and are distributed in a disordered way, 
resulting in a excellent match with experimental diffrac- 
tion data[25, 26 . Moreover, they hinder the mobility of 
Li interstitials, which would be very high if the clusters 



FIG. 2. (color online) Comparison between the intramolecu- 
lar peak of the N-H radial distribution function, as computed 
from molecular dynamics (MD) and harmonic lattice dynam- 
ics (HLD) 30 , with and without considering nuclear quan- 
tum effects. Note that because of strong anharmonicities the 
classical HLD provides an unsatisfactory description of this 
system. A similar discrepancy is found when comparing the 
quantum thermostat results and HLD with Bose-Einstein oc- 
cupations. 
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were broken [25]. Starting from this structure we per- 
formed Born-Oppenheimer molecular dynamics simula- 
tions within Density Functional Theory (DFT) with gra- 
dient corrected exchange and correlation functional [27] as 
implemented in the CPMD[28 package. Technical details 
of the calculations are the same as those used in Ref.[25 . 
Thermal averages were performed over 15 ps of trajec- 
tory, following 5 ps used for equilibration at T = 300 K. 

Nuclear quantum effects were treated by means of a 
quantum thermostat, which is based on a bespoke gen- 
eralized Langevin equation containing correlated noise. 
This stochastic process is designed to mimic the quan- 
tum mechanical phase space distribution in the harmonic 
limit. The resulting non-equilibrium dynamics samples 
a stationary distribution which is a good approximation 
of the quantum mechanical one in fairly anharmonic sys- 
tems as well. This result is achieved without requiring 
any preliminary information, except for an upper-bound 
estimate of the stiffest vibrational mode present. We used 
the set of parameters qt-50_6, which can be downloaded 
from an on-line repository [29 . We refer the reader to 
Refs.[7J[29] for further details. 

Before discussing the comparison with the experimen- 
tal proton momentum distribution, it is necessary to 
stress that the presence of anharmonic wagging modes of 
the imide groups makes a quasi-harmonic treatment in- 
appropriate. This is apparent in figure [2j where we com- 
pare the radial distribution function of the N-H group 
as computed from molecular dynamics and from the har- 
monic normal modes. The same figure also highlights 
the importance of nuclear quantum effects, which modify 
dramatically the typical fluctuations of the imide bonds. 

While hydrogen and the stretching mode in particu- 
lar exhibit the largest deviations from classical behavior, 
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lithium and nitrogen are also light nuclei and are there- 
fore subject to nuclear quantum effects, albeit to a lesser 
extent. For instance, the average kinetic temperature 
computed during the quantum-thermostatted dynamics 
deviates from the classical value, and is 415 K for Li 
atoms, 410 K for N (which is heavier but participates 
into the stiff N-H stretching mode) and 858 K for the pro- 
tons. These deviations illustrate the importance that nu- 
clear quantum effects have in determining the properties 
of Li2NH, such as the temperature at which the transi- 
tion between the low-temperature and high-temperature 
phases occurs [26 . 

Computing the proton momentum distribution from 
the quantum-thermostatted dynamics is straightforward. 
In fact, the thermostat has been designed to yield the 
correct momentum and position distribution in the har- 
monic limit, and has proven to work also in the anhar- 
monic case. Thus we only need to collect the momen- 
tum histogram to obtain the three-dimensional PMD, 
which is plotted in figure [3j The anisotropy of n(p), 
which is a purely quantum mechanical effect, reflects the 
symmetry of the local environment of protons inside the 
crystal, with the bonds aligned along (111) directions of 
the nitrogen antifluorite sublattice, pointing towards Li 
vacancies [25j [26l 02]. While we could not measure the 
directionally-resolved n(p) because of the difficulties in 
obtaining a single-crystal sample of appropriate dimen- 
sions, this result demonstrates the ease by which this 
quantity can be accessed by our computational technique, 
providing detailed information which can help to inter- 
pret experiments for other systems [32 . 

To compare with the powder sample experimental 
data, we spherically averaged the three dimensional 
PMD. As it can be seen from figure [4j there is a sat- 
isfactory agreement between experiments and theory. In 
particular, we observe that quantum-thermostatted sim- 
ulations match experimental data better than the results 
from the quantum harmonic approximation. To quan- 
tify the discrepancy with experiments, and to extract 
information relevant to the structural model of Li2NH, 
we fitted the experimental and theoretical PMD's with 
a anisotropic Gaussian model. We assumed a different 
spread in the directions parallel and perpendicular to the 
N-H bond[33ll34], resulting into 
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which was then spherically averaged. 

The resulting fit matches both curves very well, with 
values of adjusted R 2 greater than 0.999. The fit yields 
the parameters an = 6.49 A -1 and <j± = 3.15 A -1 for the 
theoretical PMD, and cry = 5.77 A -1 and a± = 3.70 A -1 
for the experimental one. 

The discrepancy is certainly larger than the experi- 
mental error bar, however the anisotropy in the PMD 



FIG. 3. Three-dimensional proton momentum distribution 
for the low temperature phase of lithium imide. Isosurfaces 
enclose 95%, 90%, 50% and 10% of the probability den- 
sity. The arrangement of the hydrogens around a Li vacancy, 
aligned along the (111) axes, is also reported, relative to the 
Cartesian reference. 
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is correctly captured, and the spread in the direction 
parallel and perpendicular to the bond is qualitatively 
reproduced. This is indeed a remarkable result for an 
approximate model of nuclear quantum effect such as 
the one used in the present work, which can be applied 
with no computational overhead with respect to stan- 
dard ab initio molecular dynamics. While it is difficult 
to assess the uncertainty in our theoretical results, a 
possible approach to gauge the error is to repeat sim- 
ulations with a different set of noise parameters. We 
did so using the parameters qt-20_6[29], and obtained 
cry = 6.31 A" 1 and a ± = 3.40 A" 1 . The compari- 
son with the results obtained using qt-20_6 sets a lower 
bound for the error at about 10%. Although the dif- 
ference between the PMD obtained for different struc- 
tures proposed for Li 2 NH|25j [26j [35j [36] is smaller than 
our error bar above on the spherically- averaged distribu- 
tion, the difference in structure reflects in qualitatively- 
different three-dimensional n(p). However, for the rea- 
sons discussed in our previous work [25 , the other pro- 
posed structures have to be dismissed because of worse 
agreement with diffraction data. 

In conclusion, we have shown how a recently-developed 
method to compute nuclear quantum effects can be used 
together with ab initio molecular dynamics to model 
complex materials containing light atoms. We were able 
to use a large simulation cell, which was mandatory for 
Li2NH to reproduce the experimental structure of the 
low-temperature phase. We obtained good agreement 
with the experimental proton momentum distribution, a 
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FIG. 4. (color online) Comparison between the spherically- 
averaged proton-momentum distribution expected for a classi- 
cal system and from harmonic lattice dynamics at T = 300 K, 
the PMD measured from a lithium imide sample and that 
computed from quant um-thermostatted molecular dynamics. 
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result of great significance given the growing importance 
of inelastic neutron scattering experiments as a sensitive 
probe of the local environment in hydrogen-containing 
materials. Both experimental and theoretical data are 
perfectly fitted by a model in which imide groups per- 
form hindered librations. 

The possibility of treating delicate nuclear quantum 
effects inexpensively, albeit approximately, suggests that 
the quantum thermostat should be used whenever light 
atoms are present, and a more accurate treatment by 
path integral methods is unfeasible because of the ex- 
cessive computational effort. Together with accurate 
first-principles calculations of the interatomic forces, this 
will shed light on the role of nuclear quantum effects in 
condensed-phase systems. R. Senesi and J. Mayers are 
gratefully acknowledged for suggestions during data re- 
duction and analysis. We thank C. Andreani for useful 
discussions. One of the authors (AP) acknowledges the 
CNISM-CNR research program. 
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